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Abstract 

To investigate the finite time singularity in three-dimensional (3D) Euler flows, the simplified 
model of 3D axisymmetric incompressible fluids (i.e., two-dimensional Boussinesq approximation 
equations) is studied numerically. The system describes a cap-like hot zone of fluid rising from 
the bottom, while the edges of the cap lag behind, forming eye-like vortices. The hot liquid is 
driven by the buoyancy and meanwhile attracted by the vortices, which leads to the singularity- 
forming mechanism in our simulation. In the previous 2D Boussinesq simulations, the symmetricial 
initial data is used, see, e.g., 20]. However, it is observed that the adoption of symmetry leads to 
coordinate singularity. Moreover, as demonstrated in this work that the locations of peak values for 
the vorticity and the temperature gradient becomes far apart as t approaches the predicted blow-up 
time. This suggests that the symmetry assumption may be unreasonable for searching solution 
blow-ups. One of the main contributions of this work is to propose an appropriate asymmetric 
initial condition, which avoids coordinate singularity and also makes the blow-up to occur much 
earlier than that given by the previously simulations. The shorter simulation time suppresses 
the development of the round-off error. On the numerical side, the pseudo-spectral method with 
filtering technique is adopted. The resolutions adopted in this study vary from 1024 2 , 2048 2 , 4096 2 
to 6144 2 . With our proposed asymmetric initial condition, it is shown that the 4096 2 and 6144 2 
runs yield convergent results when t is fairly close to the predicted blow-up time. Moreover, as 
expected the locations of peak values for the vorticity and the temperature gradient are very close 
to each other as t approaches the predicted blow-up time. 
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I. INTRODUCTION 



The formation of finite time singularities in the three-dimensional (3D) incompressible 
Euler equations is a controversial subject in fluid mechanics. There are two main reasons for 
studying the singularity development: firstly, the verification of the finite time singularity 
may aid to the understanding of the onset of turbulence in slightly viscous flows, and secondly 
if such singularities do exist then they may provide a means by which energy cascades to 
and concentrates on small scales. 

The equations under consideration are the following: 

u t + u • Vu + Vp = 0, (1) 
V-u = 0, (2) 
u(x,0) = u , (3) 

where, u is the velocity, p is the pressure, u is the given initial velocity which is smooth in 
the sense that the data reside in some Sobolov space. However, whether the solution u can 
remain smooth for all time has yet to be proved. 

Since the question of global existence of the 3D incompressible Euler equations has been 
inaccessible analytically, the question has been investigated numerically in the past two 
decades. The most commonly used quantity in determining global existence of Eqs. 
is the vorticity uj = V x u, which has been first developed by Beale, Kato & Majda [l| 
(also see jgj), and later refined by two other groups 0,0]- It is shown in |l| that u will blow 
up at a finite time T c if and only if 

uj\L°ods — > oo, as t / T c . (4) 

This result is especially useful when the question of global existence is numerically investi- 
gated, because if we find the maximum norm of the vorticity behaves like (T c — t)~ a with 
a > 1, a finite time singularity has then developed. Based on this theory, three kinds 
of numerical efforts have been made to search for the singularities in the 3D Euler flows 
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1. The original 3D Euler simulation does not adopt symmetric assumption. This scheme 
requires the largest computer resource. 



2. The symmetry introduced by Taylor & Green |5( (see also |6j for some detailed dis- 
cussion) uses only 1/64 of the total computational time for the non-symmetric case. 

3. An even further symmetric technique introduced by Kida requires only 1/192 of 
the computational time for the non-symmetric case. 

In contrast to the 3D research, 2D study is much easier to be performed analytically or 
numerically. Under the axisymmetric assumption, the 3D Euler equations can be replaced by 
2D Boussinesq convection equations. This assumption is an even more aggressive symmetric 
assumption than the Taylor-Green or Kida flows because it turns this 3D problem into a 2D 
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one, which can save computer resource significantly. We will give a brief discussion about it 
in Section II. Similar to the 3D theory, the blowup judgment has been developed for the 2D 



Boussinesq convection flows 2(], |21|, |22j . It has been realized that if the maximum absolute 
values of the vorticity and temperature gradient behave like (T c — t)~ a and (T c — t)~P with 
a > 1 and (3 > 2, a finite time singularity will be developed. This simplified 2D model, 



although lacking of the vortex reconnection as in the 3D simulations 23l. l24l. 12511 . can also 
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reveal certain possibility of singularity formation 

Among the numerical efforts in addressing the singularity issue, some highest possible 
resolutions have been attempted in the last two decades. There are mainly two kinds of 
numerical schemes involved: pseudo-spectral methods and finite differences. For the pseudo- 
spectral efforts, Kerr used the resolution up to 1024 x 256 x 128 without any symmetric 
assumption 0], and in a recent work of Hou and Li an extremely fine grid of resolu- 
tion 1536 x 1024 x 3072 is used. The Taylor-Green or Kida simulations have reached the 
effective resolution of 1024 3 12] and 2048 3 Q. Moreover, the most intensive axisymmetric 
simulation [2(J uses a resolution of 1500 2 . 

To the best of out knowledge, the most intensive 3D finite difference calculation was 
performed by Grauer et at Q. They did not make any symmetric assumption, but used the 
adaptive mesh refinement (AMR) method to enhance the effective resolution to 2048 3 . For 
the 2D Boussinesq simulation, uniform grids of 512 2 and 1792 x 1280 were used in 20j and 



[301 ] , respectively. Pumir & Siggia [2|| employed an adaptive 256 2 grid to achieve a resolution 
of 10 7 in both dimensions. A less aggressive adaptive grid, maybe more accurate due to less 
frequent re-meshing, uses a 512 2 deformed grid to reach a 4600 2 effective resolution j^ . 

Besides the numerical efforts mentioned above, the 3D Taylor series analysis, although 
limited due to the lack of proper parallelized softwares handling high precision calculations, 
is also used to analyze the 3D Euler singulari ty p roblem E T ^ 

461 ] , but the original 



analytical studies have been carried out j^Es, 3^ 4ol 
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3D Euler equations have been modified which makes the use of relevant mathematical tools 
possible, Likewise, some 2D modified Boussinesq equations were also employed in simulations 
and analysis 0, 48, 49, 5(|. Although these models do not have much connection with 



the original 3D Euler equations, they may provide some hints to the understanding of 



the singularity issue. One of the latest results was given by Frisch et al. [5jJ, |52| who 
conducted the 2D calculations in the so-called parareal domain by taking the advantages of 
both spectral and adaptive methods. 

One of the main purposes of this work is to investigate how to set up an effective initial 
condition which can be used in the simulations of the Boussinesq equations. The reason 
that we emphasize the importance of proper initial conditions is that a poorly chosen initial 
condition may not lead to blowup or may lead to blowup at a numerically unacceptable 
large time. 

Section II contains a brief discussion on our numerical method. We adopted the so-called 



phase-shifted technique to do the de-aliasing [53J, |54j. It is clear that to a reliable and 
efficient numerical is very important in simulating possible singular behaviors. Section III 
gives a quite complete discussion of three initial conditions, including subtle but significant 
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FIG. 1: The analogy between the (r, z) plane in 3D axisymmetric flows with swirl and the (x,y) 
plane of 2D Boussinesq convection flows. We keep u = near y = and y = 2tt throughout our 
2D Boussinesq doubly-periodic calculation to avoid the coordinate singularity in the corresponding 
3D axisymmetric flows (see Fig. Eta)). The solid black arrow at the bottom indicates the direction 
of the gravity, opposite to the direction of buoyancy. 

differences from earlier work. In section V, we use the parallel strategy combined with the 
traditional parallel FFT and task distribution schemes jHH to solve the 2D Boussinesq 
equations with resolutions up to 6144 2 . The corresponding effective resolution in the 3D 
Euler formulation is 7r x 6144 3 ~ 9000 3 , which is much finer than any previous efforts. 
Singularity development will be demonstrated by considering several physical quantities 
including the peak vorticity and temperature gradient. 

II. THE NUMERICAL SCHEME 

The connection between the 3D axisymmetric Euler flow with swirl and the 2D Boussinesq 
convection has been established in |2|. For completeness, we will review the connection 
briefly. 

The basic vorticity equations for the axisymmetric swirling flow in the (r, 6, z) cylindrical 
coordinates are of the form 




(5) 



(6) 



where 



D 



d r d z d 



(7) 




Dt 



dt dr dz 

dv r dv z 



(8) 



dz dr 
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It should be realized that the (z, r) plane is essentially two dimensional flows with (x, y) 
coordinate system (Fig. By introducing a streamfunction ip, the 2D inviscid Boussinesq 
equations have the following form: 

e t + u ■ ve = o, (9) 

u t + u • Vw = -9 X , (10) 
Aip = -u, (11) 

where the gravitational constant is normalized to g = (0, —1), 6 the temperature, u = (u,v) 
the velocity, u> = (0, 0, to) = V x u vorticity, and ip the stream function. It should be pointed 
out in the remaining of the work 9 is different with that indicated in it denotes the 

temperature. Comparing with the pure 2D Euler equations with the streamfunction-vorticity 
uj-ip formulation, there is an extra temperature equation, i.e. Eq. (J9~|). and an extra term X 
associated with the Buoyancy in Eq. (JTUJl . Actually, if lj° in Eqs. (JSHHJ) is replaced by u, 
and (rv e ) 2 by p, a link can be established if we evaluate all external variable coefficients in 
Eq. (0BJ) at r = 1. It should be noticed that this link leads to some coordinate singularities 
under the cylindrical coordinates. Consequently, it is essential to keep the corresponding 2D 
Boussinesq solutions away from the horizontal boundaries (here, in this paper, y = and 
y = 2n; see Figs. [Hand Fig. ED^a)). 

The method used in our numerical simulation is pseudo-spectral approximations with 



some proper de-aliasing technique. We also adopt the filterin g te chnique |59|] (a careful 
discussion on filter in turbulence simulations can be found in |60j) to modify the Fourier 
coefficients such that the stability of the numerical scheme is enhanced. The machine accu- 
racy of our computer with double precision is e = 10~ 16 ~ e -37 , and the modifying factor in 
the filter is <f(k) = e - 37 ( 2fe / 7V ) 16 for k < N/2, where N is the Fourier modes in each direction. 

We now briefly discuss the effects of different numerical schemes for the current research. 
As to be shown later, the main concern is whether or not a 5 or 5-like function can be well 
resolved by using a discrete scheme. It is well known that the spectral coefficients of the S 
function is constant for different modes: 



2 2 

6(x,y)~C J2 E e- i(nx+mv \ (12) 

N 



2 m ="T 



where C is a constant and will be set to 1 in the following. If the resolution iV 2 — > oo, 
we will get the exact Fourier representation. Of course, it is impossible to do this by any 
computer, and what we can do is to use some finest possible grid (the largest iV 2 appearing 
in the literature is limited to 8192 2 so far [57]). The de-aliasing technique is to remove the 
aliasing error by setting some of the Fourier coefficients to zero. The filter we adopted is to 
make some of the Fourier coefficients smaller. So actually, the computer-represented value 
of the S function is determined by the resolution iV 2 , the filter and the de- aliasing. 

The de-aliasing scheme to be used is the so-called phase-shift scheme 01 , which retains 
about 7/9 of the total modes. The filter narrows the gap of the peak values of de-aliasing 
pseudo-spectral schemes and pure spectral schemes. In Table I, the filtering 5 function value 
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TABLE I: Peak values when a 5(x, y) function is represented by different schemes with the same 
resolution N 2 . 





peak value before 
filtering 


percentage retained 
by the filter 


peak value after 
filtering 


phase-shift 


2vriV 2 /9 


81.3% 


0.568iV 2 


No de-aliasing 


N 2 


59.6% 


0.5967V 2 



(a) RUN A: t =3.16 



(b) RUN B:t=2.2 




(c) RUN C:t=2 
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FIG. 2: The vorticity contours of three runs when non-velocity crosses the horizontal boundary, 
a) and b) are before the blow-up time for RUN A and B, and c) is after the blow-up time of RUN 
C. 



of phase-shift scheme is only 2% lower than the no de-aliasing scheme. In our calculation, we 
follow the tradition of our pseudo-spectral code (W. H. Matthaeus, private communication) 
to use the circular truncation in our running 58]. 



III. INITIAL CONDITIONS 

In this paper, three different initial conditions will be considered, which will be named 
RUN A, RUN B and RUN C, respectively. RUN A adopted the same initial condition as 
used in Ref 0: 

w(ar,2/,0) = 0, (13) 
9(x,y,0) = 509 1 (x,y)9 2 (x,y) [1 - 9 1 (x,y)} , (14) 

where if S(x, y) := ix 2 —y 2 — (x— n) 2 is positive, 9i = exp (1 — n 2 /S(x, y)), and zero otherwise; 
if s(y) := \y — 2ti\ /1.957T is less than 1, 9 2 = exp (1 — (1 — s(y) 2 ) -1 ), and zero otherwise. By 
choosing the initial conditions (fTB^ - lfPfj) . we can test our discretization schemes by comparing 
our numerical results with those given in j2o| . 

This governing equations (P|)- (jll|) with the above initial condition have been studied 
intensively by E & Shu [2(J with spectral methods (resolutions: 1500 2 ) and ENO finite 
difference methods (resolutions: 512 2 ). They have predicted that the side part of the rising 
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(a) RUN A: t =0 



(b) RUN B:t=0 




2 3 4 5 
(c) RUN B: t =1.2 



2 3 4 5 
(d) RUN C: t =0 




FIG. 3: a), b) and d) are contour plots of temperature for three different initial conditions, d) is 
obtained after the intermediate results of RUN B (c) are compressed by the factor of 2. 

bubble is much more dangerous than the front of the bubble (Fig. Efa)). However, this 
initial condition has two shortcomings: 

• It introduces another symmetry assumption besides the axisymmetric assumption de- 
scribed in the previous section: the symmetry respect to x = ir in the (x, y) plane. 
This will cause a further deviation from the original 3D non-symmetric Euler flows. 

• It is demonstrated that the blow-up time of this flow is after t = 3.16, at which both 
y = and y = 2tc are crossed by non-zero velocity, see Fig. EJa). This is due to the 
assumption we adopted to simplify the 3D axisymmetric flow to the 2D Boussinesq 
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flow. In other words, what we are studying here can not be closely related to the 3D 
Euler flow anymore. 

To fix the first difficulty above, we change the initial temperature field in (J 14)) to 

50 

6(x, y, 0) = —(4a; - 3tt)0i(z, y)9 2 (x, y) [1 - x {x, y)\ . (15) 

TC 

The factor (Ax — 3tc) / tc introduced here is to break the symmetry assumption with respect 
to x = tc in the (x, y) plane. However, this simulation (denoted as RUN B) will also cross 
both y — and y = 2tc (see Fig. |2Jb)) before T c . Therefore, a better initial condition (RUN 
C) is to compress the intermediate results at t = 1.2 obtained in RUN B. More precisely, we 
let 

u(x,y,0)=u'(x,2y-0Aic,1.2), 9(x,y, 0) = 6'(x, 2y - 0.4tt, 1.2), (16) 

for (x, y) G [0, 2iz] x [0, tc] (where 6' and lu' are obtained by solving Eqs. (JSJ-JHJ), (|T3*jl and 
()15|) with a 2048 2 grid), and zero otherwise. The initial condition associated with RUN C is 
demonstrated in Figs. Efb, c, d). 

There are three main advantages for using the above initial conditions with RUN C: 

• There will have no symmetric and boundary crossing problems as observed in RUN A. 
The flow pattern does not cross the horizontal boundary until t = 2 (Fig. |2Jc)) while 
the predicted blow-up time is around t = 0.91 (see Section IV B). 

• Compared with RUN A, RUN C only needs about 1/4 simulation time to reach the 
blow-up time T c , and the advantages of higher resolutions show up at early stages of 
the simulations. As a result, RUN C can save quite large amount of computational 
time, which is particularly important in this kind of study. 

• Shorter simulation time also suppresses the development of round-off error, which is 
non-trivial in the current high resolution simulations (see Appendix). 

IV. NUMERICAL RESULTS 

In this section we will present a detail discussion of the numerical results of RUN A and 
RUN C. Some features of Run B can be represented by RUN A and C, and will be only 
briefly discussed. RUN A will serve as a validate run. Moreover, it is used to illustrate some 
pitfalls we should avoid, such as round-off error, symmetric effect, and coordinate singularity 
etc. RUN C is the simulation from which we will draw main conclusions. 

A. Numerical results for RUN A 

For RUN A, three resolutions are used: 1024 2 , 2048 2 and 4096 2 . The corresponding time 
steps used are 0.00004, 0.00002 and 0.00001 respectively, given by the CFL condition. 
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FIG. 4: The 3D perspective plots of vorticity for RUN A with a 4096 2 grid at different times. All 
3D effect pictures in this paper (Figs. HI ISland lllj) are based on smooth flow field. 



t = 2.0 



t = 3.0 



FIG. 5: The 3D perspective plots of temperature for RUN A with a 4096 2 grid. 



t=2 



t =2.5 



t =3.1 



t =3.2 



t=3.4 



t =3.5 



t =3.6 
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FIG. 6: (Color online). Zoom-in contour plots of vorticity in RUN A at different times with the 
resolution of 4096 2 . Only details near the predicted singularity location ([4.5,5.5] x [2, 27r]) are 
shown. The plus symbol indicates the location of |w| ma x in the whole [0, 2ir] x [0, 2tt] domain. 
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t=2 t =2.5 t =3.1 t =3.2 t =3.5 t =3.6 




3.5 4 4.5 5 3.5 4 4.5 5 3.5 4 4.5 5 3.5 4 4.5 5 3.5 4 4.5 5 3.5 4 4.5 5 

FIG. 7: (Color online). Zoom-in contour plots of temperature in RUN A at different times with 
the resolution of 4096 2 . Only details in [3.25,5.3] x [0, 2ir] are shown. The "*" symbol indicates 
the location of \V8\ max in the whole [0, 2tt] x [0, 2tt] domain. 

At the beginning of the simulations (for all three resolutions), the cap-like initial tempera- 
ture field (Fig. Efa)) will rise and develop into a "two-eye" system. Driven by the buoyancy, 
the vorticity field will also develop from the initially unified zero to a two-eye system, but 
with one positive and one negative "eyes" this time. Figs. |U and El show the 3D pictures of 
vorticity and temperature field in the whole domain, and it is clear that the peak values of 
vorticity are located along the edge of the "eyes." 

The system is basically symmetric respect to x = n if round-off error does not play a role, 
so there will have two Ic^max locations: one positive and one negative. It is worth mentioning 
that round-off error will spoil this symmetry when very high resolutions are adopted (see 
the discussion in Appendix). Our highest resolution for RUN A is 4096 2 , and the symmetry 
is well maintained before t = 3.6. After t = 3.6, the filter spoils the singularity-forming 
mechanism, and the values of |o;| maa; begin to drop when t becomes larger. In the following, 
we will only use the positive half of the vorticity field (i.e. x G [tt, 2tt]) for the divergence 
analysis. 

The "+" symbol in Figs. El indicates the locations of |o>| maa; at different times. The "+" 
is around y = 3 in the beginning, and suddenly jumps to y ~ 4.5 after t = 3.1. The flow 
field has no significant change around t = 3.1, and the front of the bubble rises continuously 
as what happens before that time. 

The "*" symbol in Figs. [7| indicates the locations of | V9\ max at different times. It should 
be noticed that the tail left behind the rising front is forming an "eye" after t = 2.5, and 
there is a smooth filament connecting the "eye" and the head of the bubble. 

The "*" symbol jumps to the edges of the forming "eye" around t ~ 3.6, roughly at the 
same time when the original smooth filament breaks up into many even smaller "eyes" (see 
Figs. ITo^) . This is because our simulation begins to become under- resolved. The filter begins 
to dramatically reduce the peak values of those 5-like functions. The "*" symbols can not 



10 



0.5 1 1.5 2 2.5 3 3.5 4 

time 
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time 
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FIG. 8: RUN A: Time evolutions of the rela- 
tive errors in T 2 , T4 and 9 max (solid line: 1024 2 ; 
dash line: 2048 2 ; circle line: 4096 2 ), defined as 
(T 2 (0) - T 2 (t))/T 2 (0), (T 4 (0) - r 4 (t))/T 4 (0) and 

\ u max 

(o) - 9 max (t)\/\e max {o)\. 
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0.5 1 1.5 2 2.5 3 3.5 4 
time 

represent the exact maximum locations of the 2D Boussinesq solution anymore. The edge 
of temperature contour is very smooth before the collapsing time, and the front part of the 
rising bubble and the two eyes left behind (Figs. EJ) stretch the filament connecting them. 

Questions arise as to how effectively the large local quantity (in our case, the vorticity and 
temperature gradients) can be resolved and how much the results vary with resolution. It 
is evident that for the extreme hypothetical case, when one of the local quantities is a delta 
function, it can only be resolved with infinite resolution. To resolve an ideal shock wave 
(with zero thickness), for example, in the absence of any smoothing, needs infinite number 
of grid points. In the following, we will check several other aspects on the effectiveness of the 
numerical simulations, mainly by finding some quantity properties useful in demonstrating 
the effect of the mesh refinement, or identifying the trend when high resolutions are adopted. 

First, we check the following three values which are time independent due to the 
divergence-free constraint, the doubly-periodic condition and the inviscid transport equation 



(Eq. ©): 
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FIG. 9: RUN 
A: Time evolu- 
tion of \uj\ m ax 
and \V0\max 

(solid line: 
1024 2 ; dash line: 
2048 2 ; circle 
line: 4096 2 ). 
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time 



FIG. 10: RUN A: The distance between the locations of |w| m ax and \V9\ max at different times for 
three kinds of resolutions (solid line: 1024 2 ; dash line: 2048 2 ; circle line: 4096 2 ). 
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• T 2 (t) = 8 2 (x,y,t)dxdy, 

• T 4 (t) = Cf^9\x,y,t)dxdy, and 

• ^ max (t) 

Fig. IS] shows that the global average quantities are well conserved within 1% error for 
4096 2 and 2048 2 resolutions throughout the entire simulation period t G [0,4]. To save the 
computation time, our 4096 2 run starts from t = 2.0 using the intermediate result of 2048 2 , so 
for 4096 2 run, those errors are defined as (T 2 (2.0) - T 2 (t))/T 2 (2.0), (T 4 (2.0) - T 4 (t))/T 4 (2.0) 
and \9 max (2.0) - e max (t)\/\9 max (2.0)\. 

The 1024 2 run has a very poor performance after t = 3.8 because none of these relative 
errors are below 10%. The performance of 2048 2 run is almost as good as on 4096 2 grid by 
considering T 2 and T 4 , with under 1% relative errors. However, the relative errors on 9 max of 
4096 2 are always below 10~ 4 , which are only about 2% of the errors in 2048 2 simulation after 
t = 3.5 (Figs. Et c ))- These errors give some indicators on how far the numerical solution is 
away from the real one. 

Figs. IHJa) and (b) show that T 2 and T 4 errors of 4096 2 run are always lower than those 
given by the 1024 2 and 2048 2 simulations. However, due to the round-off error in our 
calculations (see Appendix), this is not the case for 6 max error (Figs. Efc)), because the 
errors associated with the 4096 2 resolution seems to be at the same level of those with the 
2048 2 resolution from t = 2.0 to t = 2.6. When t is close to the blow-up time, the 4096 2 
run conserves the value of 8 max much better than other runs. In fact, when t — > T c , the 
truncation error will overwhelm the influence of machine precision, and higher resolutions 
have more advantages than lower ones. For future simulations, however, we may have to 
take some measures to control the round-off error when refined grids are used, say, the 
quadruple precision (the machine accuracy e = 10~ 32 ) may be used to replace the present 
double-precision (e = 10~ 16 ). On the other hand, the present results of 2048 2 and 4096 2 
simulations seem to be accurate enough to draw some conclusions in the blow-up analysis. 

Figs. Efa) and (c) show the time evolutions for |u;| maa ; and |V#| maa; with different resolu- 
tions. We re-plot these maximum value evolutions in Figs. 0(b) and (d) with logarithmic 
scales. It seems that the growth of | oj\ max and \W9\ ma x in the 4096 2 run terminates in a 
finite time with I 

cu\max ~ ('/*c — i) and crudely |V^| m aa; ~ ('/*<: — t) ■ The correspond- 
ing predicted blow-up time is T c = 3.72. Even for 2048 2 run, there is an obvious trend 
towards the singularity. 

However, there are still two facts that cause suspects of this conclusion: 

• The simulation never reaches the blow up time T c ; 

• Even for the 4096 2 run, we can not say that we have resolved the problem confidently 
because the time-evolution curves for | climax and iV^m^ show quite large differences 
between high and low resolutions at late stages of the simulations. 

The first point is unavoidable due to the discrete nature of the numerical simulations. The 
second point may be improved by employing higher resolution, or, at least we can extend 
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the overlapping region to a time closer to T c . For example, the \uj\ max (t) lines of 1024 2 and 
2048 2 runs diverge around t = 2.7, but the 2048 2 and 4096 2 curves diverge around t = 3.2. 
However, the effort in this direction may be useless because, as explained in the following 
paragraph, we can not get a clearer physical picture of the singularity forming mechanism 
even we use higher resolution. 

Physically speaking, if a singularity is about to form at t = T c on the point (x c ,y c ), 
the locations of |co>| max and |V0| maa ; should approach (x c ,y c ) when t — > T c . This means the 
distance between the locations of these two peak values should become smaller and smaller 
as t — > T c . Although it is impossible for any numerical solver to reach T c , we can use this 
property to check whether a numerical simulation reflects the physical mechanism properly. 
Fig. ITU1 shows that the locations of |o;| max and | V8\ ma x in RUN A initially stay far away from 
each other (about 1.7 length unit in the [0, 2n] x [0, 2n] domain) and they do not get closer 
when the "singular" time is approaching 67]. The distance are noticeably increased after 
t « 1.7 for all three resolutions. When t G [3.2,3.5], it is even wider (see also the specific 
locations of those maximum values in Figs. Eland[7j). The physical mechanism behind it is 
very unclear, maybe because the numerical code is trying to preserve the symmetry from 
the very beginning of flow pattern (Chapter 2 of t 66]), accompanied by the singularity 
forming mechanism. It seems that there are at least two locations at which the singularity 
is formed simultaneously at x G [ir, 27r], fighting for the locations of global maximum values. 
Because the divergent conclusion drawn from Figs. Mjo) and (d) means nothing unless higher 
resolutions are adopted, we should consider some other initial condition which may reveal 
more physical phenomena more effectively. The results of RUN A seem too complicated to 
be analyzed. 

To sum up for this subsection, we basically end up with similar conclusion as in [2^ 
in the sense that the edge of the "cap" is more dangerous than the front. Although our 
code can resolve the 2D Boussinesq equations accurately until about t « 3.5 on the 4096 2 
grid, we can not obtain a clear physical picture relevant to the singularity issue for the 3D 
Euler equations. RUN A can only be used to reveal some singular phenomena in the 2D 
Boussinesq equations. 

B. Numerical results for RUN C 

In RUN C, we use four sets of grids, namely 1024 2 , 2048 2 , 4096 2 and 6144 2 . The corre- 
sponding time steps are 0.00004, 0.00002, 0.00001 and 0.000008 respectively (for 4096 2 run, 
we use At = 0.00002 for t < 0.22, and 0.00001 for t > 0.22 to save the computation time). 

Fig. ^2 shows the 3D pictures of vorticity in the whole domain of RUN C. The rising 
bubble also develops into two "eyes" like that in RUN A. In RUN C, as a whole, the positive 
"eye" has much larger absolute values for |u;| and |V0| than those of the negative "eye." 
There is no symmetry with respect to x = it anymore. In the following, we will only exam 
the details around the positive "eye" with a particular attention to the locations of the 
global maximal |a>| and |V0| (see Fig. fl2"j) . 

A noticeable difference between Fig. Upland Fig. |H]is about the location for |o;| maa; . For 
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FIG. 11: The 3D perspective plots of vorticity for RUN C with the 6144 2 grid. 
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FIG. 12: (Color online). Zoom-in contour plots of vorticity (the first row) and temperature (the 
second row) in RUN C at different times with the resolution of 4096 2 . The "+" and "*" indicate 
the location of |to>| maa; and |V#| maa; | in the whole [0, 2tt] x [0,2-/r] domain, respectively. 

RUN C, from the very beginning until the collapsing time (t ~ 0.88), the location of |co> (man is 
changing continuously and smoothly without any noticeable jumping. This is quite different 
to the sudden jumping of "+" symbol in RUN A at t ~ 3.2. The \W6\ max is located around 
x ~ 4 in the beginning (until t ~ 0.16, see "*" in Fig. H2|L right on the head of the rising 
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bubble (in this sense, it behaves like RUN A). 

The "*" jumps to x ~ 5 at t ~ 0.18, roughly at the same time when the "eye-like" vortex 
begins to form. Afterwards, the "*" symbol moves quite continuously towards the location 
of |o;| max . After t = 0.7, the locations of |^| ma x an d |Vt9| maa; are very close to each other. 
In the meanwhile, the filament connecting the cap front and the "eye" in the temperature 
field becomes thinner and thinner, which breaks down after t = 0.87 due to insufficient 
resolutions. 

For the three low resolutions (1024 2 , 2048 2 and 4096 2 ), the time evolutions of the three 
time independent values T 2 (£), T 4 (t) and 6 max {t) in RUN C (Fig. H3|) are quite similar to the 
results in RUN A (Figs. [SJ). Because shorter simulation time is needed, the 4096 2 result of 
RUN C is less affected by the round-off error. As a result, the 4096 2 curve is always below 
those of two lower resolutions, see Fig. ITHlf c). The truncation error is the main source of 
the error for these three low resolutions. 

Furthermore, if we we only look at errors of T 2 and T 4 , the 6144 2 results show much 
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better conservation property than that of the three coarser grids (Figs. IT3T a) and (b)). 
The 9 ma x error, which is more easily affected by the round-off error, shows no advantage of 
finer grids (Figs. H^T c)): for < t < 0.6 the error of the 6144 2 resolutions is at the same 
level as that of the 4096 2 resolution. However, when t > 0.6, the advantage of the 6144 2 
resolution is observed: the error of 9 max remains around 10 -6 which is one order lower than 
that of the 4096 2 result. The round-off error of 6144 2 can roughly match the truncation 
error. This indicates that it seems necessary to employ quadruple precision instead of the 
double precision if extremely finer resolutions (say 10 5 in each dimension) are used. 

Again, we investigate the time evolution of the distance between the locations of |o;| max 
and | V6* | max under the same resolution. Unlike the RUN A case (see Fig. [TUj) . we obtain 
a very clear physical picture (Fig. [T4^) . Around t = 0.2, the distance experiences a sudden 
drop, and then decays slowly. The 1024 2 curve has an early dramatic increase at t m 0.65 
due to under-resolution. The results of the 4096 2 and 6144 2 resolutions are quite satisfactory 
since the distance remains small. The finest resolution shows no obvious advantage here. For 
RUN A, the first-time jump of "+" (Figs. EJ) is followed by the increasing distance between 
the locations of |cG>| ma:2 ; and \W9\ max (see Fig. ITU)) . For RUN C, however, it seems that the 
singularity forming mechanism is the main driving source, which leads to the decreasing 
distance between the locations for l^max and \VQ\ max after the first jump of "+". 

Now we can analyze the singularity forming mechanism of RUN C. We focus on the 
filament connecting the cap front and the right "eye-like" vortex (see, for example, the 
temperature contour plot at t = 0.8 in Figs. [T2"j) . The filament corresponds to the hotter 
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region of the fluid field, which will rise due to the buoyancy without any other forces. On 
the other hand, the forming "eye" tries to absorb the fluid around it. These two mechanisms 
fight with each other, which makes the filament thinner and thinner until singularity appears. 
For RUN A, besides these two mechanisms, there is a third one joining the competition, 
namely, the symmetry with respect to x = it. The third mechanism is not physical, but the 
numerical scheme tends to preserve it throughout the simulations. With three mechanisms 
working simultaneously, it is difficult to know when and where the singularity will be formed, 
although the time evolutions of |co>| ma:r and |V#| ma a: give some hints on solution blow-up. In 
contrast, the locations of |co>| ma 2; and |V0| ma a: in RUN C are getting closer and closer when 
the solutions become singular. 

We now think it is safe to carry out some singularity analysis for RUN C. Figs. fToT a) 
and (c) show the time evolutions for |co>| maa; and |V#| maa: in different resolutions; the inverse 
of the maximum values are plotted in Figs. IToT bVdnd ). In the following, we will pay our 
attention to the 4096 2 and 6144 2 curves, because these two lines seem to be overlapping 
until around t = 0.88. 

For the 4096 2 run, there are some very small vorticity structures which begin to appear 
at t = 0.84 in the lower part of the smooth outer layer where the maximum \u\ turns up. 
For 6144 2 run, this happens at t > 0.87. After this critical time, the filter we adopted in the 
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codes removes more and more energy from the system. Consequently, the global average 
values like T2 are greatly affected. As a result, our numerical results after t > 0.87 may not 
be reliable. Actually, it is observed that the |a>| moa! and | V9\ max experience a drop-down after 
t ~ 0.9. This indicates that the filter prevents the occurrence of the blowup of the maximum 
vorticity, which is similar to the viscosity effect in high Reynolds number simulations (e.g. 



12j). For the 1024 2 and 2048 2 runs, however, the drop-downs appear later than t ~ 0.88 
and the blow-up time T c is also delayed. 

Only the sample maximum values before t = 0.86 in 6144 2 run will be used in the 
singularity analysis because before this critical time convergence between the 4096 2 and 
6144 2 results is observed. With a statistical weighted least square fitting for the data after 
t = 0.6, we conclude that the growth of |o;| max and | V6\ max for RUN C terminates in a finite 
time with I Climax ~ (T c — j cLIld. CFud-Gly ~ (T c - t)~ z ' iS with T c = 0.91. 

There seems to have a trend of singularity also for the 2048 2 run, and the 1024 2 result 
indicates no blow-up. Moreover, the collapsing times for higher resolutions are earlier than 
those of lower ones because the filter removes more energy with low resolutions, retarding 
the singularity forming mechanism. 

It should be noticed that the peak vorticity we obtained is the value modified directly 
by the filter, and the highest peak of the 5-like function (see the peak near "+" in Fig. fTTI 
when t = 0.86) is reduced significantly. On the other hand, |V#| max is less affected by the 
filter because it is the temperature field being filtered (not the V# field). No place in the 
temperature field is really close to a 5-function. Hence, the temperature filed is less affected 
by the filter than the vorticity field. From this point of view, the |V#| max curve is more 
accurate than the |co>| ma 2; one. If there is a blow-up in the simulation, the |V6*| maa: curve will 
show a stronger divergent tendency than the |k>| maa ; one (Figs. fTo"]) . 

We have not performed the simulation on grids finer than 6144 2 , but we can predict 
some results from the present computations. It is well known that when a delta function 
is approximated by a finite number of Fourier modes, each doubling of the resolutions will 
cause doubling of the maximum value and 2 n+1 times the maximum values of the n th -order 
space derivative (see the analysis in Appendix A of 12j). In the final state of our simulation 
(Fig. H2](t=0.86)), the cut-line at y = 2.72 through the out-layer of the "eye" of the vorticity 
field looks very similar to a delta function. Therefore, when finer and finer resolutions are 
used the peak vorticity and temperature gradient are getting larger and larger. From this 
point of view, further 2D Boussinesq simulations with even higher resolutions will support 
our singularity prediction although it seems impossible for any code to reach the same T c if 
the current filtering relevant scheme is used. 



V. DISCUSSIONS AND CONCLUSIONS 

There have been extensive discussions on the 3D Euler singularity with viscous simula- 
tions by using the Taylor-Green vortex and high symmetry flows (the most intensive one, in 
our opinion is ^^)- As the Reynolds number is increased, the amplitudes of the maximum 
vorticity, skewness, and flatness increase, and the peaks are attained at earlier times. This is 
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quite similar to what happens to our simulation with the filter. Our simulation also reaches 
an earlier "blow-up" when higher resolutions are adopted. These tendencies give hints that 
for the pure inviscid case (no viscosity or filter), all these quantities may blow up. 

A natural question is to discuss what is the influence with symmetry assumptions. We 
tend to believe that for the Taylor-Green, high-symmetry Kida and axisymmetric flow, all 
the symmetry constraints are unstable, and in the absence of symmetries, the flow will 
escape from the singularity formation direction. Moreover, singularity formations depend 
strongly on how the initial conditions are set up. The numerical method chosen has also a 
strong impact on the results. For example, finite difference simulations that do not make use 
of any symmetry produce no sing ularity indication , whereas Fourier- Chebyshev !()■ 
simulations with symmetry constraints sug gest singular trends. 

To sum up, we follow the track of jij and end up with a more developed result. The 
disagreement between the locations of |co>| m ax and | V6\ max suggests that the symmetric initial 
data proposed in 2(| may not be appropriate (see Fig. fTUj) . To fix this problem, a new initial 
condition is proposed in this work, which enables us to make some fine grid simulations. 

There are several points that make us believe that there is a singularity in RUN C: 

• Our simulations show that the time evolution curves of |w| ma x and \VQ\max become 
steeper when the grids are refined. 

• The distance between |u;| max and |V#| ma;c on the 6144 2 and 4096 2 grid shows a con- 
vergent behavior near T c . 

• The distance between the locations of |^| ma x and |V#| max is small even up to the 
predicated blow-up time. 
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APPENDIX: THE EFFECT OF FILTER ON ROUND-OFF ERROR 

The influence of round-off errors on computation has been observed and noticed for a long 
time, probably starting from the appearance of the digital computer. The round-off error 
starts to occur at the very beginning of all numerical simulations, and can be amplified by 
the time integration procedure. A comprehensive study on how machine precision can affect 
the dynamic simulation is carried out in [61j . Later, in vortex sheet roll- up simulations, a 
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practical result about these effects with two kinds of precisions (7 and 14-digit arithmetic) 
is presented j^. In the following, we will try to give a brief analysis about the round-off 
error in our simulations. 

There are three ways to accumulate round-off errors in our simulations: 

1. Fast Fourier Transform (FFT) can control the initial round off error very well. Through 
the result in j6^,|6^, we can roughly guess that the error is amplified about 100 times 
for current resolutions (1024 2 — 4096 2 ), that is, the error will be in order of 10 -13 if 
the double precision (e = 10~ 15 ) is adopted. The 4096 2 grid will have only slightly 
larger round-off error than that of 1024 2 grid (only a factor of /og4096/Zogl024 m 1.2, 
based on the worst possibility given in 63]). 

2. The filtering technique we adopted can significantly amplify the round-off error. Ac- 



cording to 65[, the filter will amplify the round-off error by a factor of (J2 c|) 2 . Here, 
Ck is the factor added onto the Fourier coefficients by the filter. Actually, most of Ck 
are 1, so roughly speaking, the round-off error will be enlarged by N times for a A^ 2 
simulation. This means our 1024 2 run will have a 10 -10 error for each time step, and 
the error in 4096 2 runs is four times larger. 

3. Remember that the analysis above only happens within each time step, and the real 
round-off error can be even larger in a dynamic run. The time step of 4096 2 run is 
1/4 of the 1024 2 run according to the CFL condition, which means 4096 2 needs four 
times as many time steps as 1024 2 . 

As a whole, for the same simulation, the round-off error in 4096 2 run will be 1.2 x 4 x 4 = 
19.2 times larger than that in 1024 2 run. This difference is already big enough to break the 
initial symmetry maintained by the numerical solver. Fig. IToT a) is the contour plot of RUN 
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A at t = 4.0 with the 1024 2 grid, which has less details and level than the 4096 2 grid (Fig. 
ITffl b)). However, Fig. IT6T a) preserves the symmetry respect to x = ir very well, while Fig. 
[Tffl b) has a different number of vortices on the left and right sides of the picture. In the 
main text, a discussion of the round-off error is also carried out together with Figs. IHland 
IT3~l focusing on the round-off error amplified by the time evolution. 

Figs. EH are results after the simulation becomes under resolved. Our 4096 2 simulation 
of RUN A is very symmetric respect to x = n before t = 3.5. The round-off error in 
6144 2 will be about 30 times larger than that of the 1024 2 run. The round-off error is not 
negligible comparing with truncation error for 6144 2 run (Fig. ITBT c)). it seems a must to 
adopt quadruple precision on a finer grid. 
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